function [t, y] = rk4(f, t0, y0, tf, dt)  
    % RK4 算法求解常微分方程  
    % 输入参数：  
    % f    - 函数句柄，表示 dy/dt = f(t, y)  
    % t0   - 初始时间  
    % y0   - 初始状态向量  
    % tf   - 最终时间  
    % dt   - 时间步长  
    % 输出参数：  
    % t    - 时间向量  
    % y    - 状态向量（每一时间步长的解）  
  
    % 初始化时间向量和状态向量  
    t = t0:dt:tf;  
    y = zeros(length(t), length(y0));  
    y(1,:) = y0;  
  
    % RK4 算法迭代计算  
    for i = 1:length(t)-1  
        t_current = t(i);  
        y_current = y(i,:)';  
  
        % 计算 k1, k2, k3, k4  
        k1 = f(t_current, y_current);  
        k2 = f(t_current + 0.5 * dt, y_current + 0.5 * dt * k1);  
        k3 = f(t_current + 0.5 * dt, y_current + 0.5 * dt * k2);  
        k4 = f(t_current + dt, y_current + dt * k3);  
  
        % 更新状态向量  
        y(i+1,:) = y_current + (dt / 6) * (k1 + 2 * k2 + 2 * k3 + k4);  
    end  
end 